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In this paper we study the interplay between stochastic gene expression and system design using 
simple stochastic models of auto-activation and auto-inhibition. Using the Poisson Representation, 
a technique whose usefulness in the context of non-linear gene regulation models we elucidate, we 
are able to write down exact results for these feedback models in the steady state. We exploit this 
representation to analyze the parameter-spaces and demarcate where different behaviors including 
power-law, conventional bimodal, a novel bimodal with graded characteristics and sub-Poisson noise 
occur. Using our results, we reexamine how well the auto-inhibition and auto-activation models serve 
their conventional roles as paradigms for noise suppression and noise exploitation respectively. 
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I. INTRODUCTION 

Stochastic fluctuations occur inevitably in the bio- 
chemical reactions that implement cell functions result- 
ing in cell-to-cell variability of the components, in partic- 
ular, in protein numbers. This noise in protein numbers 
can either be functionally significant or may be detrimen- 
tal and thus require suppression [1]. Such fluctuations 
arise even in a population of cells that were initially iden- 
tical due to the inherently probabilistic nature of the key 
biochemical processes involved in gene activation, tran- 
scription and translation coupled with the small numbers 
of molecules involved [2, 3]. Stochastic gene expression 
and speciflcally, such 'intrinsic' fluctuations in protein 
numbers, have been the focus of several experimental and 
theoretical studies [2-10]. An important systems biolog- 
ical quest is to understand if the activity of a protein is 
determined predominantly by just its abundance (aver- 
age value) or if its cell-to-cell variability (noise strength 
and shape) is also an important determinant, given a 
gene network architecture. 

Previous studies have indicated that tight control of 
protein numbers, when desired, is often achieved by an 
auto-repression motif of gene expression [11-15]. It has 
even been argued that the reason why this motif oc- 
curs far more frequently in nature (40% of the known 
transcription factors in E. Coli are controlled by neg- 
ative auto-regulation [16]) than in studies of random- 
ized networks is because it achieves stability against fluc- 
tuations [17]. On the other hand, noise may be ex- 
ploited by cells to switch between different expression 
states especially in key cellular processes where 'lock- 
ing' of sub-populations of cells into distinct fates needs 
to be achieved without changing the underlying network 
structure. The auto-activation motif has been implicated 
in such systems when population heterogeneity is desir- 
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able [18-21]. 

It is now well established that in many systems pro- 
teins are produced in 'bursts' with significantly varying 
dead times in between bursts. This important aspect of 
gene expression is encapsulated in a model in which the 
gene can stochastically switch between long-lived "off" 
states and "on" states leading to spurts of mRNA and 
protein expression. The dynamics of protein bursting de- 
pends on the interplay of time scales, the switching rate 
and degradation rates, in such a model. We explore the 
interplay between stochastic gene expression and system 
design by studying two simple stochastic bursting models 
with auto-regulation, i.e., with feedback from the protein 
produced to the gene switching rate. In the positive (neg- 
ative) feedback model, the amount of protein produced 
proportionally increases the propensity of the gene to be 
in the "on (off)" states. 

We use the Poisson representation, a technique whose 
usefulness in the context of solving nonlinear models of 
gene expression we elucidate, to derive not just the exact 
steady state protein distributions in these models, but 
to also analyze the parameter spaces and to demarcate 
where in them bimodal, power-law tailed, sub-Poisson 
and other distributions occur. Using these results we 
re-examine how well these two models, characterized by 
negative and positive feedback, serve as paradigms for 
noise suppression and noise exploitation respectively. 

While the idea of writing down protein distributions 
as exact linear superpositions of Poisson distributions is 
relatively new [9] , mixtures of Poisson distributions have 
long been studied in various contexts including photon 
statistics in quantum optics [22] and in accident prone- 
ncss models in actuarial sciences [23]. Remarkably, in 
both the auto-activation and the auto-repression cases, 
we find classes of mixed Poisson distributions that, to the 
best of our knowledge, have not been previously consid- 
ered [24] . What is more, they arise dynamically in these 
models. We also show that the Beta-Poisson mixture, 
which has been prc;viously studied versatile prior 
distribution in accident proneness models [23], arises as 
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a limiting case from these distributions. 



II. THEORETICAL FRAMEWORK 

Detailed expositions of the Poisson Representation can 
be found in [23, 25] amongst others. We have briefly dis- 
cussed the Poisson Representation for linear models of 
gene regulation without feedback such as the gene puls- 
ing model, in [9]. We analyze how the exact steady-state 
protein distributions P{p) of the models with feedback 
discussed in this paper may be represented as a superpo- 
sition of Poisson distributions, with a weighting probabil- 
ity density p(A) for the Poisson mean A. In other words, 
wc examine whether a probability density p(A) can be 
found such that 
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If such a p{X) is foimd, an immediate consequence is 
that the corresponding P{p) is super-Poisson, i.e., has 
a Fano-factor (FF) > 1. p{X) is a function of a contin- 
uous variable in contrast to P{p) and its convexity and 
monotonicity properties are easier to ascertain and these 
in turn determine the shapes of P{p). Bimodal P{p) dis- 
tributions correspond to concave(upwards) p(A); Power- 
law tails in P{p) arise when p(A) itself has a monotoni- 
cally decreasing power law region, and a monotonically 
increasing p(A) leads to unimodal P(p) distributions with 
a mode around the upper edge of the A interval. When 
p(A) is concave downwards with a maximum at some in- 
termediate value of A, then unimodal P(m) distributions 
with a mode around the same value obtains. We use the 
exact, analytical expression that wc derive for p{X) to 
map out where in the parameter space each ciualitatively 
distinct shape arises. 

For the auto-activation model schematically shown in 
Figure SI, let PQ{p,t) and Pi{p,t) denote the probabili- 
ties that there are p proteins at time t and that the gene 
is in the off and on state respectively. Master Equations 
describing the time-evolution of these probabilities can 
be written down using a standard procedure [25]. 

Let us define po{X) and pi{X) as 
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Note that p{X) — po{X)+pi{X) satisfies the normalization 
condition / dXp{X) = 1. The Master equations for A) 
are give by 

dtPo{X,t) = -Cfpo{X,t) + CbPi{X,t) + dx[Xpo{X,t)] 
- a (Apo - dx{Xpo)) (3a) 

dtPi{X,t) = CfPo{X,t) - Cbpi{X,t) + dx[Xpi{X,t)] 

+ a {Xpa - dx{Xpo)) - PbdxPi{X,t) , (3b) 
with the boundary condition 
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FIG. 1: Top Left: The response of the protein distribution to 
increasing activation strength in the (j) < 1, (5 < 1 quadrant 
in the auto-activation model resembles the classic 'binary re- 
sponse' associated with auto-activation systems. Bottom Left: 
The effect of increasing activation strength, a, on the pro- 
tein distribution in the fourth quadrant, i.e., 4> < > 1 
is 'graded' in this quadrant. Top Right: Effect of increasing 
auto-repression strength, r, on the Fano-factor of the protein 
distribution in the auto-repression model. For the values cho- 
sen, both /5 and the Fano-factor go through maximum values 
at (different) intermediate values of r, before the distribu- 
tion becomes sub-Poisson after the threshold value r = ro; 
this happens exactly when ^ = 0. Bottom Right: Effect of 
increasing auto-repression strength r, on the protein distribu- 
tion in the auto-repression model. Six different points from 
the above figure arc chosen from the range where fS remains 
positive. The other rates are the same as the ones used in the 
previous figure, r increases from lighter to darker values. 



for < A < Amax) where a Xmax needs to be computed. 
If we can solve the above equations we have p(A) and the 
Poisson representation exists by construction. 

Similarly, for the auto-repression model shown in Fig- 
ure SI, if the Poisson representation exists, we can find 
the corresponding densities Pa(A) that satisfy the Master 
Equations 

dtPo{X,t) = -Cfpo{X,t) + CbPi{X,t) + dx[Xpo{X,t)] 
+ r (Api - <9a(Api)) (5a) 

dtPi{X,t) = CfPo{X,t) - CbPi{X,t) + dx[Xpi{X,t)] 

- r (Api - dx{Xpi)) - pbdxPi{X,t), (5b) 



with the boundary condition 
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for < A < Xmax, where such a A^ 



must be found. 
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III. RESULTS 

To place our results for the auto-activation and auto- 

rcprcssion models in context, wc will find it useful to 
compare these results with those derived for the linear 
pulsing model in [9] . Both auto- regulation models reduce 
to the linear pulsing model (LPM) in the limit where 
the auto- activation strength, a, or the auto-repression 
strength, r, tends to 0. In [9] we have shown how the 
'phase-diagram' of all possible distributions for the linear 
pulsing model can be drawn in terms of the two scaled 
rates Cf/pd and Cb/pd- 

A. Auto-activation 

For the Auto- activation model (Figure SI) we have 

found a solution to the Master Equations (3), showing 
the validity of the superposition-of-Poissons representa- 
tion, given below: 

p{X) = X^~' - xY'*" (7) 

with < A < Pb/pd', A/" is the normalization con- 
stant. This exact cxprc;ssion k;a(ls naturally to the cor- 
rect parametrization of the combinations of the rate 
constants that are relevant for analyzing this nonlinear 
model. We rescale A by Pb/pd so that it lies between 
and 1. It is useful to rescale all rates by the effective 
protein degradation rate, Pd- For convenience in clas- 
sifying the different kinds of protein distributions that 
arise in this model, we define the following parameters: 
a = apb/{l + a), (j) = c//(l -|- a) and /3 = Cb/{1 + a). We 
then have 

p(A) =Afe"^ A'^-i (1 - Xf-^ (8) 

Note that cj) and /3 characterize the singularity at the 
upper and lower limits of A. Using this superposition- 
of-Poissons representation we have found that in each of 
the four 'quadrants' determined by (f> and /3 greater or 
less than 1, the protein distribution has a distinct shape. 
Since the superposing density, p{X), is found to extend 
from X = to X = pb, P{p) extends till ^ pb- When the 
density diverges at both limits, i.e., (p and /? < 1 yielding 
a p that is concave upwards the protein distribution is 
bimodal. When p vanishes at both limits, i.e., ^ and 
/3 > 1, yielding a p that is concave downwards a broad 
bell-shaped distribution of proteins arises. 

As the autoactivation strength a— >-0, a— ^0, (^— ^c/ 
and /3 — >■ Cf, the protein distribution of the autoactivation 
model becomes the exact steady-state distribution [9] ob- 
tained in the LPM. The latter is a Beta distribution 

p{X) =Af X^i-^ {l-Xy"-'^ (9) 

Thus, the 'phase-diagram' of possible distributions in this 
model is very similar, in large regions of the parameter 



space cf) and /3, to that of the LPM, despite the auto- 
activation, once we identify <p and /3 in this model with 
Cf and Cb in the simple pulsing model. 

Wc focus on the most interesting new feature that 
arises in this 'phase-diagram' in this model. Consider 
the quadrant where cj) < 1 and /3 > 1. When a = 0, i.e., 
in the LPM, we have found [9] long-tailed distributions 
with power-law behavior. In the auto-activation model, 
in contrast, two possibilities arise depending on whether 
a is lesser or greater than ac = (Vl — <P+ VP ~ 1 
the former case long-tailed distributions with power-law 
regions arise, with an exponent — l asm the a = case. 

For a > Uc the distribution becomes an unusually be- 
haved bimodal distribution! To appreciate its nature we 
recall that when both cf) and 13 are < 1 (Figure 1) bi- 
modal distributions occur with the two modes always at 
and Pb, i.e., at the edges of the allowed values of A. As 
a the activation strength increases, the weights aroimd 
and Pb , are redistributed without affecting the separation 
between the modes. This is the classic 'binary response' 
[18] typically associated with auto- activation: cells may 
be thought to be divided into two sub-populations with 
low and high protein numbers and increasing activation 
strength only changes their relative proportions. In con- 
trast, the new bimodal distribution exhibits a second 
mode not at pb, the maximum allowed value of A, but 
at intermediate values. As a is increased by increasing 
a, the protein distribution goes from being monotonically 
decreasing power-law to bimodal because aiito-activation 
affects cells with intermediate numbers of proteins the 
most. Thus, when the feedback strength is strong enough 
that a > ttc, a new minimum and as well as a new maxi- 
mum develop in p(A), at intermediate values of A. Corre- 
spondingly, P{p) becomes bimodal with the second mode 
arising at a value oi p < pb- As the activation strength 
increases, this mode tends to higher values of p, but the 
weight at (the first mode) simultaneously erodes rapidly 
making the distribution effectively unimodal for strong 
enough activation. Thus in this quadrant even though 
bimodal distributions arise for intermediate activation 
strength, the response to increasing activation is really 
'graded' as illustrated in Figure 1. As a increases, the 
protein distribution goes from being negatively skewed, 
with a large likelihood of obtaining a small number of 
proteins to a positively skewed distribution, with a large 
likelihood of obtaining a large number of proteins. 

We point out the possible relevance of our results to 
the observation in a recent experiment of Maheshri et 
al. [26] of bimodal protein expression in a synthetic yeast 
system with positive feedback and no cooperativity as in 
our model. As the activation increases, their distribu- 
tion goes from a broad bell-shaped distribution to the 
bimodal distribution similar to the one described above. 
Our model explains their observation of graded response 
of the IxtetO promoter with increasing auto activation 
strength. As expected from our model, with increasing o, 
the mode at larger value travels further towards the right 
and acquires more weight until a Poisson like distribution 
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occurs. 

Since this model is nonlinear the equations for all the 
moments are coupled and one needs the full distribu- 
tion to obtain even the lowest two moments. Explicit 
expressions are displayed in the Supplementary section. 
These yield the Fano-factor (FF), the ratio of the vari- 
ance to the mean of a distribution. FF may or may not go 
through a maximum value as a is increased, but beyond 
a threshold FF always decreases with a and tends to 1 as 
a ^ 00. Thus increasing auto-activation results in noise 
reduction, a role not conventionally associated with pos- 
itive feedback. This is true since the gene is always "on" 
as the activation strength tends to infinity, and a Pois- 
son protein distribution, with FF = 1 results. For any 
initial choice of parameters, for large enough a, both (p 
and j3 fall below 1, and the protein distribution becomes 
bimodal. However, in the limit a ^ oo, the mode at 
is entirely eroded and p{X) 5{ph — \): Pip) becomes 
Poisson. 



B. Auto-repression 

The analysis of this model proceeds along the same 
lines as the auto-activation model. Once again, this 
formulation leads naturally to the correct parametriza- 
tion of the combinations of the rate constants that are 
relevant for analyzing this nonlinear model. We define 
the new parameters, a = rpb/{l + r)^, (f) = Cf and 
/3 = pbr/{l + rY + (cfc — c/r)/(l + r). All rates have 
been scaled by the protein degradation rate, pd as be- 
fore. In terms of these new variables, the steady state 
generating function is identical in form to that derived in 
the auto-activation model! 

However, there is a subtle difference which has pro- 
found consequences: the parameter (3 can become neg- 
ative for suitably chosen rates, pb,Cf,Cb and r, in this 
model, unlike the auto-activation case. Thus the weight- 
ing probability density p{X) can be found, only if /3 > 0. 
This immediately implies that for ^ > 0, the protein dis- 
tribution in the auto-repression model is super-Poisson, 
i.e., its FF is > 1 and thus "noisier' than the Poisson 
distribution that arises in the simple birth-death model. 

When /3 is < 0, we find that the protein distribution 
becomes sub-Poisson, i.e., its FFs becomes < 1. Thus, 
only when /3 < can the auto-repression be said to be 
strong enough to cause reduction of the noise level in re- 
lated models, such as the LPM and the auto-activation 
model. On analyzing the condition ^ < 0, we find that 
for any given value of the rates c/,Cb and pb, there is a 
threshold value of the repression strength, ro, such that 
when r increases beyond this threshold value, the distri- 
bution becomes sub-Poisson (as illustrated in Figure 1). 
As seen in the top right panel of Figure 1 this suppres- 
sion occurs for large values of r and over a narrow range. 
Exactly at the threshold value, the Fano-factor is found 



to be unity. The expression for ro is 

^ Cb+Pb-Cf + ^/{cb+Pb - CfY -I- 4c/ Cj, 

2c/ ^ ' 

For values of r > ro, i.e., when the distribution is sub- 
Poisson, a formal expression for p{\) may be derived, 
with the understanding that it can no longer be inter- 
preted as a probability density. In fact, A now extends 
over the complex plane. Remarkably, even in this case, 
the functional form of p(A) remains the same for a suit- 
ably chosen contour. In this case, depending on whether 
C/ is < 1 or > 1, the protein distribution is a monoton- 
ically decreasing or a sharply peaked bell-shaped distri- 
bution, respectively. 

When /3 > 0, we find that the different possible distri- 
butions of the auto- activation models all occur for auto- 
repression for appropriate values of a, {3 and v, when 
^ > 0, as illustrated in Figure 1. This underscores the 
inadvisability of naively inferring that the choice of auto- 
inhibition motif is designed to obtain noise suppression 
without further exploring the specific details of the sys- 
tem. See also [27] for a control and information theoret- 
ical perspective on the issue. Quantitatively, A is in the 
range 0<A<pb/(l-|-r) and so the protein distribution 
extends to about p ~ pf,/(l -t-r). The efl'ective parameter 
/3 is now a function of all the rates in the problem while 
the effective parameter (/> = c/ as in the linear pulsing 
model. 



IV. CONCLUDING REMARKS 

The auto-regulation motif is ubiquitous in gene regu- 
lation [3, 17]. The auto-regulation models studied here 
are admittedly simplified descriptions of those observed 
in nature: we have not included separate transcription 
and translation steps. In prokaryotes, since mRNAs are 
rapidly translated into proteins, this is typically a reason- 
able approximation. For eukaryotic systems, when the 
mRNA time-scale is significant, these results should not 
be applied literally. Further, the effects of co-operative 
auto-regulation are not included in our models. However, 
even in this simple model a plethora of behaviors are ob- 
served including power laws, and bimodal distributions 
that behave in a graded fashion and sub-Poisson statis- 
tics. We have also established the utility of the Poisson 
Representation which yields quite naturally, the impor- 
tant, scaled, dimensionless parameters that characterize 
non-linear gene regulation models. We saw that auto- 
activation produces 'binary' responses to increasing acti- 
vation strength and that auto-repression produces noise- 
suppressed sub-Poisson protein distributions in very lim- 
ited regions of the parameter space. Our work serves to 
add a note of caution to assuming that positive and neg- 
ative feedback, when found in natural biological systems, 
are present to serve these purposes. 
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Supplementary Section for 
Mixed Poisson distributions in exact solutions of 
Stochastic Auto-regulation Models 
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o4 1 Derivation of steady state distributions 

^ 1.1 Introduction 

^1 We introduce the method that we use to investigate the auto-activation and auto- 

Q repression models in the context of the hnear pulsing model(LPM) where there is no 

• ^ feedback. This model is described by the reactions 

o 

oo 

p\j Let the probabilities that the gene is in the on or off state at time t and exactly 
p proteins are present be denoted by Pi{p,t) and Po{p,t), respectively. The Master 

^ Equations satisfied by Pq and Pi can be shown to be the following [1] . 

Tl dPoip,t) 



pHo (1) 



.'^ dt 

>< dPi{p,t) 



B dt 



-CfPo{p,t) + ChPi{p,t) + ka[{p + l)Po{p + l,t)-pPo{p,t)[ 
CfPo{p,t) - CbPl{p,t) + kd[{p+l)Pi{p+l,t)-pPi{p,t) 

+h[plip-l,t) - Plip,t)\ 



One can solve this model using generating functions [1]. They are defined by G{z, t) = 
^2"^^^ P{p,t) where P{p,t) = Po{p,t) + Pi{p,t). Wc solve the models with feedback 
using the Poisson representation and we will illustrate the method in the linear model. 
We assume the existence of probability densities po{X, t) and pi(A, t) that yield Po{p, t) 
and Pi(p,t) by 

/•oo y 

Pa{p, t) = dXpa{X, t) — where a = or 1. (2) 
Jo P- 
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Figure 1: The auto-regulation models: The auto- activation/auto-repression model is 
shown on the top/bottom. D and D* denote the gene in the 'off' and 'on' states, P 
denotes the protein. In both cases, auto-regulation is assumed to use monomers of 
the protein. 



Note that the normalization condition Yl'^=o ^(P^^) ~ ^ leads to the corresponding 
condition J dXp{\) = 1 on p{X,t) = po{X,t) + pi{X,t). From the Master Equations 
for PQ{p,t) and Pi{p,t) we obtain the equations satisfied by the Poisson densities po 
and pi, 

dtpo{X,t) = -CfPo{X,t) + Cbpi{X,t) + kddx[Xpo{X,t)] (3) 
dtpi{X,t) = Cfpo{X,t) - Cbpi{X,t) + kddx[Xpi{X,t)] - kbdxPi{X,t) . (4) 

In deducing these we impose the condition that the boundary terms resulting from 
integration by parts vanish: 



kd {Xmax - A) —7 Pi (A) 
ml 



(5) 



where A lies in the interval < A < Xmax- The boundary terms vanish identically at 
Xmax, provided that Pi(A) is less singular than l/iXmax — X) and at the lower limit if pi 
vanishes at A = 0. The solution we obtain pi(A) does indeed behave as required and 
so the assumption that the boundary terms vanish can be justified a posteriori. We 
also have Xmax — physically reasonable since when the DNA is always on it 

leads to a Poisson distribution for the protein distribution with the A-parameter equal 
to the maximum rate of protein production kb/k^. In the steady state, the coupled 
Master Equations for po and pi can be solved and yield 



p(A) = X-^+'f (kb - X)^-^+'"'^ 



(6) 
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where the normahzation factor M can be easily found by integrating over the allowed 
range of A. Since pi cx \p the vanishing of the boundary terms is easily justified. The 
result obtained is that in steady state the distribution p(A) is a scaled Beta Distribu- 
tion with Beta parameters c/ and c;, and with A ranging from to /cf,, i.e., Xmax = ^fe- 
The steady state distribution P(p), give by a superposition of Poisson distributions 
with the probability p(A) obtained above; using the integral representation of the 
confluent hypergeometric function 



^0 



we obtain 

P{p) oc iFi(c/+p,c/ + C6+p;-/c6), (8) 

where the proportionality constant can once again be found from normalization. We 
follow this procedure to determine without recourse to the standard generating 
function method for the models with positive and negative feedback below. 



1.2 The Auto Activation Model 

This model is described by the following reactions with the protein switching the gene 
from the off- to the on-state: 

D%D* (9) 

Cb 

D + P ^D* + P (10) 
D* J^D* + P (11) 

(12) 

The protein probability distributions with the DNA in the D or D* states denoted by 
Po{p,t) and Pi{p,t) respectively satisfy the Master equations 

^^^^^ = -CfPoip,t) + c,Piip,t) - apPoip,t) 

+ k, [{p + l)Po(p + 1, t) - pPoip, t) ] (13) 

dPi{p,t) 



dt 



CfPo{p,t) - ChPi{p,t) + apPo{p,t) + ka[{p+l)Pi{p + l,t)-pPi{p,t)] 
+ kt,[P,{p-l,t) - P^{p,t)] (14) 
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Figure 2: Steady state p(A) distributions for the LPM; they are labeled by (cfe,c/) 
in units of fc^. The translation rate is lOOfc^ for all the distributions. The figure 
shows prototype distributions for the four major regimes Cf,Cb ^ fc^ in parameter 
space corresponding to the four quadrants. The case where all the regimes overlap 
with the uniform distribution for Cf = Cb = is also shown. All distributions extend 
from A = to A = fcfe = 100. 



As before, we can rewrite the Master Equations for P„(p, t) in terms of po and pi. 
The result is 



dtpo{X,t) = -c/po(A,t) + CbPi(A,t) + kddx[Xpo{X,t)] 

- a (Apo - dx{\po)) 
dtpii\,t) = c/po(A,t) - CbPi(A,t) + kddx[\pii\,t)] 

+ a (Apo - <9a(Apo)) - hdxpi{\t) . 



with the requirement that the boundary terms vanish, i.e. 

^max 



± a e '^^ A po(A, t) 
m\ 



A'^ 







(15) 
(16) 

(17) 



where A is restricted to the interval < A < \max- Once again, the solution we find 
will satisfy these boundary conditions at A = and Xmax = k^/kd- Expressing all 
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rates in units of kd we can solve the steady-state Master Equations to obtain 

p(X) ^ ^^e^^X(^''^ (h- X)(^~^ (18) 

with A e [0, kb]. We can use the probabihty density of p(A) to deduce the steady-state 
P{p) directly without the use of generating functions, as before. We find 

P{P) oc iF, (^+P, ^ +P; -j^) , (19) 
\l + a 1 + a 1 + aj 

where the proportionality constant is determined from normalization. 

1.3 The Auto Repression Model 

This model is described by the following reactions: 

D^D* (20) 

Cb 



D* + P—^D + P (21) 

D* D* + P (22) 

P . (23) 
We can derive the Master Equations satisfied by the A-densities as before. 

dtPo{X,t) = -CfPo{X,t) + CbPi{X,t) + kddx[Xpo{X,t)] 

+ r (Api - dxiXpi)) (24) 

dtPi{X,t) = c/po(A,i) - CbPi{X,t) + kddx[Xpi{X,t)] 

- r (Api - dxiXp,)) - kbdxPi{X,t) . (25) 



In this case, the requirement that boundary terms vanish translates to 



{kb - kdX - rX) — pi{X) 
ml 



(26) 



If we set Xmax = kb/{l+r) (using rate constants measured in units of kd) the boundary 
terms can be made to vanish provided pi(A) is well behaved at the limits. On solving 
these equations, we find that 

pi(A) oc e^r^ X^f {kb- X{l + r)f-^ (27) 

where we have defined 

o _ Cb-rcf + kbr/{l + r) 

^= YT~r ■ ^^^^ 
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For /3 > 0, the solution is well behaved and the boundary conditions are satisfied; thus 
the description of the protein distribution as a superposition of Poisson representation 
is well-defined in this case. The density p(A) is found to be 



p{X) ^ J\f e^r>^ X^f-' - Xj for/3>0 (29) 

If /3 < the boundary conditions are not satisfied and the representation as a super- 
position of Poisson distributions docs not exist. Indeed, as will explicitly demonstrate 
later, P{p) becomes sub-Poisson for /3 < with a Fano factor less than 1. Therefore, it 
is not surprising that the protein distribution cannot be described by a superposition 
of Poisson distributions. However, formally we can write down a complex representa- 
tion as used in quantum optics for squeezed states. [2] This follows from the fact that 
the integral representation of the confluent hypergeometric function used previously 
to obtain P(m), 

/ dxx^-^iu-x^-^d^'' = B{fi,iy)u''+''-\Fi{u,fi + iy;f3fi) (30) 
Jo 

can be extended for < and u < when the integral is taken along the Pochhammer 
contour. Using the Pochhammer contour, p(A) retains the same functional form even 
when P{p) becomes sub-Poisson, i.e., when /? < 0, but evidently it can no longer 
be interpreted as a probability density in this case. When convergence of boundary 
terms is imposed we can obtain p(A) from pi and thence derive the steady-state P{p) 
directly as before. We find that, for all /3 



r{cf+p)r{cf + /3) f h Y 1 iFi(cf + m, p + Cf+p; 



2 Exact expressions for the mean and the variance 

Since both the auto-activation and auto-repression models are non-linear feedback 
models, it is not possible to derive the exact mean of the protein distribution from the 
rate equations, since such a description ignores the correlation between the promoter 
and protein dynamics. As is well-known because of the non-linearity, the Master 
Equation cannot be used to derive an evolution equation for just the mean (in contrast 
to linear models): the equation for the mean is coupled to higher moments leading 
to a hierarchy of coupled equations that cannot be truncated if one desires an exact 
expression even for the first moment; the entire probability distribution needs to be 
explicitly evaluated to obtain even the mean. From the exact distribution that we 
have derived we can obtain the mean and the variance. 
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2.1 The Auto activation Model 

The steady-state mean, /i, is seen to be the first moment of p(A) and is given below: 

rp f Cf ^ _ a+Cj,+ej- + l _ gfcj, \ 

„^(p)= ' ''-'V ' '-^'^ /A (32) 

We can also obtain from the second moment of p(A) 

i^i(;3T + 2;^ + 2;^)c,(a + c, + l)^,^ 

{P[p - 1)) = J r . (33) 

1^1 (c. + c,)(a + c, + c, + l) 

The variance cr^ = ijp') — fi^ can be written down from the above expressions. As 
before, we have expressed all rates in units of the decay rate ka- In the limit a — > 0, 
the model reduces to the LPM and indeed these expressions reduce to the simple 
results for the pulsing model[l], 

— — and 

Cfe + Cf 

V (c6 + c/)(cb + c/ + l); 

The first term the variance is the noise in the simple Poisson process leading to a 
Fano factor of unity while the second term in the expression above corresponds to the 
additional noise from the DNA flipping between the on and off states. 



2.2 The Auto Repression Model 

In this model, the expression for the mean is 

^' = ^ — 7 -( V\ \^ (3^) 

y' ^ ^> 1-^1 I "^Z' (r+l)2 ' (r-+l)2 ) 
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where iFi is the regularized confluent hypergeometric function. The expression for 
the second moment is 

k,c, ik, + r+ 1) L + 1; + 1; (^) 

U — ^- — 

klc, (r (c, + h- Cf) + c, + r- (-c,)) ,h (c, + 1; --^-zMir^'^^ + 2; 

(r + ij li^l I C/, '(^TTF) 

(36) 

and from these expressions, the variance cr^ can be written down using a"^ = u — fi^. 
Once again, the LPM resuhs for the mean and variance are recovered in the hmit 
when the feedback strength becomes 0, i.e. as r — )■ 0. 

3 The phase diagram of allowed steady state dis- 
tributions 

Prom the functional form of p(A), the important dimensionless combinations (ratios) 
of the rates in the model can be deduced and they can be used to demarcate the 
regions of the parameter space of the model where qualitatively distinct distributions 
occur. 

3.1 The linear pulsing model 

Once again, it will be useful to compare the results for the auto-activation and auto- 
repression models with that obtained for the linear pulsing model and so we recall 
them here. Recall that for this case p(A) = AfX-^+^f {kb - Xy-^+^'l First, the shape 
of the distribution is evidently determined by cj and Cb, is just a scaling variable 
that determines the extent of the distribution. By examining the possible extrema of 
p(A) for the shape determining variables Cf and c^, it can be seen that the parameter 
space can be divided into the four 'quadrants' shown in Figure( 2). In the first 
quadrant (top, right), p(A) has one maximum at < A < /cb and thus the resultant 
Pip) is a broad bell-shaped curve. In the second quadrant (top, left), p(A) diverges at 
X — kh and so the corresponding P{p) is a narrow Poisson like distribution with a bell- 
shaped background. Indeed, as Cb — )■ 0, the distribution becomes Poisson. In the third 
quadrant (bottom, left), p(A) is U shaped diverging at both ends with a minimum at 
< X < kb- As a result, P{p) is bimodal with modes around and kb- Finally, in 
the fourth quadrant, i.e., the bottom right, p(A) is monotonically decreasing with a 
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maximum at and a power law tail with an exponent Cf/kd — l and the resultant P{p) 
is also monotonically decreasing with a power law tail with an exponent Cf/ka — l- 

3.2 The Auto activation Model 

For this case, it was shown that the Poisson density is given by p(A) = N e (1+")^ A (1+'") 

A) and < A < /cft. The extent of the distribution is thus determined by /c^. 

By examining the extrema of the function, it can be shown that the shapes of the 
distributions in this model are determined by = c//(l + a) and /3 = £5/(1 + a) 
and a = akb/{l + a). Qualitatively, the parameter space can still be divided up into 
four quadrants in terms of {4>,/3) similar to it was in terms of (c/,Cfe) in the LPM 
case. This is shown in Figure(3). Qualitatively, the behaviors of p(A) and P{p) in the 
first, second and third quadrants remain the same as that for the LPM case, if the 
parameter space is characterized by (0, (3) in the place of (c/, Cb) in the linear pulsing 
model. In the fourth quadrant, the behavior of p(A) becomes a dependent as shown 
in Figure(4). p(A) is monotonically decreasing as before with a power law exponent 
equal to (p/ka - 1 for a < (v^T^ + y/^^f. For a = (^T^ + ^/J^f, p(A) 
develops a point of inflection at intermediate A. For a > (Vl — + y//3 — 1^, p(A) 
has two modes, one at and the other at a value of A < kb- In this case, the system 
develops a new mode at a finite value of A which is less than the edge of the interval at 
kb, in contrast to the the bimodal in the third quadrant. As the value of a is increased, 
the second mode tends toward kb- When a — )■ 00, the weight at the mode at tends 
to and the second mode tends to kb and the distribution tends to a Poisson with 
mean kb/kd- We emphasize that this new kind of bimodal is not present in the LPM 
and is a result of the feedback included in the model. ^ 
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We note that in the third quadrant, when c/, Cb are themselves < kd, true bimodal distributions 
are obtained for P{p) but if c/, C5 > kd such that (j), (3 are still in the third quadrant, then p(A) has 
a minimum very close to and so P{p) has a small delta function mode at and is nominally a 
bimodal, but for all practical purposes has a mode just at kb- 
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Figure 3: The phase diagram of the superposing density p(A) for the Auto-activation 
model, in terms of the variables and /3 defined in the text. 

[Notice how strongly reminiscent this phase diagram is of the corresponding figure 
for the LPM, in terms of the variables c/ and Ch. Since this figure is meant to be 
illustrative, the actual values of the parameters used to generate the distribution in 
each quadrant are not important. The variable A has been rescaled by PblVd so that 
the distribution extends from A = to A = 1.] 
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Figure 4: Appearance of a new kind of bimodal distribution in the Auto-activation 
model. 

[While Figure 3 shows how strongly reminiscent distributions obtainable in the Auto- 
activation model are of those in the LPM, a closer look at the 'fourth-quadrant' in 
Figure 3 reveals that a new kind of bimodal distribution appears when appropriate 
conditions are satisfied by the parameters ( a > (i/l — -f- a//? — 1 )^), while still in 
this quadrant. The leftmost distribution is the same as the one shown in the 'fourth- 
quadrant' of Figure 3. As a is increased (from left to right), such that the < 1 
and /3 > 1, a new mode develops at intermediate values of A (< 1) when the above 
inequality is satisfied.] 



